I'm looking into doing a delta_sigma emulator. This is testing if the cat side works. Then i'll make an emulator for it.



In [1]:

    
from pearce.mocks import cat_dict
import numpy as np
from os import path
from astropy.io import fits



In [2]:

    
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()



In [3]:

    
z_bins = np.array([0.15, 0.3, 0.45, 0.6, 0.75, 0.9])
zbin=1



In [4]:

    
a = 0.81120
z = 1.0/a - 1.0

Load up a snapshot at a redshift near the center of this bin.



In [5]:

    
print z









    



0.232741617357



In [6]:

    
cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[a]}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!

cat.load_catalog(a, particles=True)



In [7]:

    
cat.load_model(a, 'redMagic')



In [16]:

    
params = cat.model.param_dict.copy()
#params['mean_occupation_centrals_assembias_param1'] = 0.0
#params['mean_occupation_satellites_assembias_param1'] = 0.0
params['logMmin'] = 13.4
params['sigma_logM'] = 0.1
params['f_c'] = 1.0
params['alpha'] = 1.0
params['logM1'] = 14.0
params['logM0'] = 12.0

print params









    



{'logM1': 14.0, 'logMmin': 13.4, 'f_c': 1.0, 'logM0': 12.0, 'sigma_logM': 0.1, 'alpha': 1.0}



In [17]:

    
cat.populate(params)



In [18]:

    
nd_cat =  cat.calc_analytic_nd()
print nd_cat









    



0.000600656711146



In [19]:

    
rp_bins = np.logspace(-1.1, 1.5, 16) #binning used in buzzard mocks
rpoints = (rp_bins[1:]+rp_bins[:-1])/2



In [20]:

    
ds =  cat.calc_ds(rp_bins)



In [21]:

    
plt.plot(rpoints, ds)
plt.loglog();

Use my code's wrapper for halotools' xi calculator. Full source code can be found here.



In [ ]:

    
xi = cat.calc_xi(r_bins, do_jackknife=False)

Interpolate with a Gaussian process. May want to do something else "at scale", but this is quick for now.



In [23]:

    
import george
from george.kernels import ExpSquaredKernel
kernel = ExpSquaredKernel(0.05)
gp = george.GP(kernel)
gp.compute(np.log10(rpoints))



In [24]:

    
print xi









    



[  2.60657131e+03   1.23027666e+03   4.58736076e+02   1.49469349e+02
   4.17659588e+01   1.27924782e+01   6.01294431e+00   4.45841995e+00
   2.54847662e+00   1.36017653e+00   7.30093914e-01   3.86348464e-01
   1.74698080e-01   6.42636931e-02   2.12185559e-02]



In [25]:

    
xi[xi<=0] = 1e-2 #ack



In [26]:

    
from scipy.stats import linregress
m,b,_,_,_ = linregress(np.log10(rpoints), np.log10(xi))



In [27]:

    
plt.plot(rpoints, (2.22353827e+03)*(rpoints**(-1.88359)))
#plt.plot(rpoints, b2*(rpoints**m2))

plt.scatter(rpoints, xi)
plt.loglog();



In [28]:

    
plt.plot(np.log10(rpoints), b+(np.log10(rpoints)*m))
#plt.plot(np.log10(rpoints), b2+(np.log10(rpoints)*m2))
#plt.plot(np.log10(rpoints), 90+(np.log10(rpoints)*(-2)))

plt.scatter(np.log10(rpoints), np.log10(xi) )
#plt.loglog();









    Out[28]:





<matplotlib.collections.PathCollection at 0x7fc0e1f1a450>



In [29]:

    
print m,b









    



-2.34633655775 2.2132896



In [30]:

    
rpoints_dense = np.logspace(-0.5, 2, 500)

plt.scatter(rpoints, xi)
plt.plot(rpoints_dense, np.power(10, gp.predict(np.log10(xi), np.log10(rpoints_dense))[0]))
plt.loglog();

bias = cat.calc_bias(r_bins) print np.sqrt(bias)

# plt.plot(rpoints, bias) plt.xscale('log') plt.ylabel(r'$b^2$') plt.xlabel('r [Mpc]') plt.title('Bias, "Updated" Values') plt.ylim([0,8]);

This plot looks bad on large scales. I will need to implement a linear bias model for larger scales; however I believe this is not the cause of this issue. The overly large correlation function at large scales if anything should increase w(theta).

This plot shows the regimes of concern. The black lines show the value of r for u=0 in the below integral for each theta bin. The red lines show the maximum value of r for the integral I'm performing.

Perform the below integral in each theta bin:

$$ w(\theta) = W \int_0^\infty du \xi \left(r = \sqrt{u^2 + \bar{x}^2(z)\theta^2} \right) $$

Where $\bar{x}$ is the median comoving distance to z.



In [36]:

    
#a subset of the data from above. I've verified it's correct, but we can look again. 
wt_redmagic = np.loadtxt('/u/ki/swmclau2/Git/pearce/bin/mcmc/buzzard2_wt_%d%d.npy'%(zbin,zbin))

The below plot shows the problem. There appears to be a constant multiplicative offset between the redmagic calculation and the one we just performed. The plot below it shows their ratio. It is near-constant, but there is some small radial trend. Whether or not it is significant is tough to say.



In [37]:

    
from scipy.special import gamma
def wt_analytic(m,b,t,x):
    return W*b*np.sqrt(np.pi)*(t*x)**(1 + m)*(gamma(-(1./2) - m/2.)/(2*gamma(-(m/2.))) )



In [38]:

    
plt.plot(tpoints, wt, label = 'My Calculation')
plt.plot(tpoints_rm, wt_redmagic, label = 'Buzzard Mock')
#plt.plot(tpoints_rm, W.to("1/Mpc").value*mathematica_calc, label = 'Mathematica Calc')
plt.plot(tpoints, wt_analytic(m,10**b, np.radians(tpoints), x),label = 'Mathematica Calc' )

plt.ylabel(r'$w(\theta)$')
plt.xlabel(r'$\theta \mathrm{[degrees]}$')
plt.loglog();
plt.legend(loc='best')









    Out[38]:





<matplotlib.legend.Legend at 0x7fc0e15a3490>



In [ ]:

    
wt_redmagic/(W.to("1/Mpc").value*mathematica_calc)



In [ ]:

    
import cPickle as pickle
with open('/u/ki/jderose/ki23/bigbrother-addgals/bbout/buzzard-flock/buzzard-0/buzzard0_lb1050_xigg_ministry.pkl') as f:
    xi_rm = pickle.load(f)



In [ ]:

    
xi_rm.metrics[0].xi.shape



In [ ]:

    
xi_rm.metrics[0].mbins



In [ ]:

    
xi_rm.metrics[0].cbins



In [ ]:

    
#plt.plot(np.log10(rpoints), b2+(np.log10(rpoints)*m2))
#plt.plot(np.log10(rpoints), 90+(np.log10(rpoints)*(-2)))

plt.scatter(rpoints, xi)
for i in xrange(3):
    for j in xrange(3):
        plt.plot(xi_rm.metrics[0].rbins[:-1], xi_rm.metrics[0].xi[:,i,j,0])
plt.loglog();



In [ ]:

    
plt.subplot(211)
plt.plot(tpoints_rm, wt_redmagic/wt)
plt.xscale('log')
#plt.ylim([0,10])
plt.subplot(212)
plt.plot(tpoints_rm, wt_redmagic/wt)
plt.xscale('log')
plt.ylim([2.0,4])



In [ ]:

    
xi_rm.metrics[0].xi.shape



In [ ]:

    
xi_rm.metrics[0].rbins #Mpc/h

The below cell calculates the integrals jointly instead of separately. It doesn't change the results significantly, but is quite slow. I've disabled it for that reason.



In [ ]:

    
x = cat.cosmology.comoving_distance(z)*a
#ubins = np.linspace(10**-6, 10**2.0, 1001)
ubins = np.logspace(-6, 2.0, 51)
ubc = (ubins[1:]+ubins[:-1])/2.0

#NLL
def liklihood(params, wt_redmagic,x, tpoints):
    #print _params
    #prior  = np.array([ PRIORS[pname][0] < v < PRIORS[pname][1] for v,pname in zip(_params, param_names)])
    #print param_names
    #print prior
    #if not np.all(prior):
    #    return 1e9
    #params = {p:v for p,v in zip(param_names, _params)}
    #cat.populate(params)
    #nd_cat =  cat.calc_analytic_nd(parmas)
    #wt = np.zeros_like(tpoints_rm[:-5])
    
    #xi = cat.calc_xi(r_bins, do_jackknife=False)
    #m,b,_,_,_ = linregress(np.log10(rpoints), np.log10(xi))
    
    #if np.any(xi < 0):
    #    return 1e9
    #kernel = ExpSquaredKernel(0.05)
    #gp = george.GP(kernel)
    #gp.compute(np.log10(rpoints))
    
    #for bin_no, t_med in enumerate(np.radians(tpoints_rm[:-5])):
    #    int_xi = 0
    #    for ubin_no, _u in enumerate(ubc):
    #        _du = ubins[ubin_no+1]-ubins[ubin_no]
    #        u = _u*unit.Mpc*a
    #        du = _du*unit.Mpc*a
            #print np.sqrt(u**2+(x*t_med)**2)
    #        r = np.sqrt((u**2+(x*t_med)**2))#*cat.h#not sure about the h
            #if r > unit.Mpc*10**1.7: #ignore large scales. In the full implementation this will be a transition to a bias model. 
            #    int_xi+=du*0
            #else:
                # the GP predicts in log, so i predict in log and re-exponate
            #    int_xi+=du*(np.power(10, \
            #            gp.predict(np.log10(xi), np.log10(r.value), mean_only=True)[0]))
    #        int_xi+=du*(10**b)*(r.to("Mpc").value**m)

        #print (((int_xi*W))/wt_redmagic[0]).to("m/m")
        #break
    #    wt[bin_no] = int_xi*W.to("1/Mpc")
        
    wt = wt_analytic(params[0],params[1], tpoints, x.to("Mpc").value)   
    chi2 = np.sum(((wt - wt_redmagic[:-5])**2)/(1e-3*wt_redmagic[:-5]) )
    
    #chi2=0
    #print nd_cat
    #print wt
    #chi2+= ((nd_cat-nd_mock.value)**2)/(1e-6)
    
    #mf = cat.calc_mf()
    #HOD = cat.calc_hod()
    #mass_bin_range = (9,16)
    #mass_bin_size = 0.01
    #mass_bins = np.logspace(mass_bin_range[0], mass_bin_range[1], int( (mass_bin_range[1]-mass_bin_range[0])/mass_bin_size )+1 )

    #mean_host_mass = np.sum([mass_bin_size*mf[i]*HOD[i]*(mass_bins[i]+mass_bins[i+1])/2 for i in xrange(len(mass_bins)-1)])/\
    #                    np.sum([mass_bin_size*mf[i]*HOD[i] for i in xrange(len(mass_bins)-1)])
        
    #chi2+=((13.35-np.log10(mean_host_mass))**2)/(0.2)
    print chi2
    return chi2 #nll



In [ ]:

    
print nd_mock
print wt_redmagic[:-5]



In [ ]:

    
import scipy.optimize as op

args = ([p for p in params],wt_redmagic, nd_mock) PRIORS = {'f_c': (0, 0.45), 'alpha': (0.6, 1.4), 'logMmin':(10.9,13.6), 'logM1': (13., 14.1), 'logM0': (9,16), 'sigma_logM': (0.01, 0.5)}

results = op.minimize(liklihood, np.array([v for v in params.values()]) ,args,method = 'L-BFGS-B', bounds = [PRIORS[p] for p in params])



In [ ]:

    
results = op.minimize(liklihood, np.array([-2.2, 10**1.7]),(wt_redmagic,x, tpoints_rm[:-5]))



In [ ]:

    
results



In [ ]:

    
#plt.plot(tpoints_rm, wt, label = 'My Calculation')
plt.plot(tpoints_rm, wt_redmagic, label = 'Buzzard Mock')
plt.plot(tpoints_rm, wt_analytic(-1.88359, 2.22353827e+03,tpoints_rm, x.to("Mpc").value), label = 'Mathematica Calc')

plt.ylabel(r'$w(\theta)$')
plt.xlabel(r'$\theta \mathrm{[degrees]}$')
plt.loglog();
plt.legend(loc='best')



In [ ]:

    
plt.plot(np.log10(rpoints), np.log10(2.22353827e+03)+(np.log10(rpoints)*(-1.88)))
plt.scatter(np.log10(rpoints), np.log10(xi) )



In [ ]:

    
np.array([v for v in params.values()])

#Try integrating over z and u jointly, explicitly nz_zspec = hdulist[8] #N = 0#np.zeros((5,)) N_total = np.sum([row[2+zbin] for row in nz_zspec.data]) dNdzs = [] zs = [] W = 0 wt2 = np.zeros_like(tpoints_rm) ubins = np.linspace(10**-6, 10**2.0, 1001) for bin_no, t_med in enumerate(np.radians(tpoints_rm)): print bin_no int_xi = 0 for row in nz_zspec.data: N = row[2+zbin] dN = N*1.0/N_total dz = row[2] - row[0] dNdz = dN/dz H = cat.cosmology.H(row[1]) x = cat.cosmology.comoving_distance(row[1]) for ubin_no, _u in enumerate(ubins[:-1]): _du = ubins[ubin_no+1]-ubins[ubin_no] u = _u*unit.Mpc du = _du*unit.Mpc r = a*np.sqrt((u**2+(x*t_med)**2).value)#*cat.h#not sure about the h #print r if r <= 10**1.7: int_xi+=du*(np.power(10, \ gp.predict(np.log10(xi), np.log10(r), mean_only=True)[0]))*dNdz*dN*H*2.0/const.c wt2[bin_no] = int_xi

plt.plot(tpoints_rm, wt2, label = 'My Calculation') plt.plot(tpoints_rm, wt_redmagic, label = 'Buzzard Mock') plt.ylabel(r'$w(\theta)$') plt.xlabel(r'$\theta \mathrm{[degrees]}$') plt.loglog(); plt.legend(loc='best')

wt_redmagic/wt2



In [ ]: